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Abstract 


The aim of this project was to investigate and improve upon existing methods of 
analysing data from COMPTEL on the Gamma Ray Observatory for neutrons emitted 
during solar flares. In particular, a strategy for placing confidence intervals on neutron 
energy distributions due to uncertainties on the response matrix has been developed 
We have also been able to demonstrate the superior performance of one of a range of 
possible statistical regularisation strategies. A method of generating likely models of 
neutron energy distributions has also been developed as a tool to this end. The project 
involved so vmg an inverse problem with noise being added to the data in various ways 

To achieve this pre-existing C code was used to run Fortran subroutines which performed 
statistical regularisation on the data. H 
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1 Background 

A solar flare is a huge release of energy within the sun which results in a sudden brightening 
of a small part of the sun’s upper chromosphere. The flare’s brightness increases dramatically 
over a short period of time, usually of the order of a few minutes. The flare’s brightness then 
slowly decays over a longer timescale, usually of the order of an hour. The energy is released 
by the flare in various forms including both particles and penetrating radiation. This project 
dealt with data from neutrons emitted during such a flare. These in turn are a product of 
reactions between mildly relativistic ions and ambient ions. Thus they may reveal aspects of 
the energy and time-dependence of flare particle acceleration. 

2 COMPTEL 

OOMPTEL is the imaging COMPton TELescope currently in orbit on the Gamma Ray Obser- 
vatory. It is used to detect spectra of neutrons emitted during solar flares which have energies 
between approximately 20 and 120 MeV. 

COMPTEL uses elastic scattering to measure the incident neutron’s direction and energy 
and consists of two arrays of detectors, one at the front and one at the rear, both of which are 
surrounded by charged particle detectors. The forward detectors consist of a liquid organic 
scintillator which has a low carbon to hydrogen ratio. This is necessary to ensure that as 
many of the incident neutrons as possible scatter from hydrogen rather than carbon. These 
detectors are also designed so that the incident neutron should scatter only once before leaving 
the detector. The telescope can also be used to detect 7-rays but it is relatively simple to 
discriminate between neutrons and 7-rays because neutrons scatter protons in the first detec- 
tor and 7-rays scatter electrons. These give rise to current pulses of different durations in the 
charged particle detectors. It is also possible to measure the time-of-flight of the neutron or 
photon from the forward detectors to the rearward detectors which also helps to distinguish 
between them as neutrons will take a longer time to reach the second detectors than photons 
travelling at the speed of light. The rearward detectors consist of Nal and are also surrounded 
by charged particle detectors. An ideal telescope event is when the neutron scatters elastically 
on hydrogen only once in a forward detector and is fully absorbed in a rearward detector. 

The amount of energy deposited by the neutron in the forward detector E\ can als o be 
obtained from the current pulse generated in the charged particle detectors. The energy, E 2 , 
that the neutron retains after it has been scattered in the first detector can be inferred from 
the time-of-flight between forward and rearward detectors. By adding the two amounts it is 
possible to calculate the total energy of the incident neutron. 1 

Etotal — E\ 4 - E 2 

If we assume that the neutron is scattered elastically, then from kinematics we can deduct 
that the scatter angle of the neutron is given by 

4> = tan ~ i (E l /E 2 ) 

where <p is the scatter angle. 2 This allows us to check that the incident neutron has been 
emitted by the sun because if we check the position of the sun at the moment when the 
neutron was recorded, it ought to appear in or near an imaginary cone whose sides lie at angle 
<j>. See Figure 2. 

All of the above is true if the neutron scatters elastically from hydrogen. Of course, other 
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Figure 1: Checking the sun’s position 
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events are possible, for example, if the neutron scatters from carbon rather than hydrogen, but 
these are recognisable from the shapes of the current pulses in the charged particle detectors 
which are very nearly 100% efficient at identifying charged cosmic rays. This project assumed 
that the data being analysed had been “cleaned up” and only results which had been identified 
as elastic scattering events were included. However, it would be possible to use the methods 
described in this report to analyse data which included obviously erroneous results. However 
a different choice of statistical regularisation method may have to be made 

COMPTEL is unique in that it has a high signal-to-noise ratio for neutron detection, 
giving it greater spectroscopic capabilities than have previously been available. It also has a 
field-of-view of 90° for neutron detection, which is quite large compared to previous detectors. 
For solar flares this should increase the number of flares which give usable data. It also opens 
up possibilities for studying Earth ‘albedo’ neutrons. 


3 Statistical Regularisation 

The task of analysing data received at COMPTEL in order to reproduce what was actually 
emitted at the sun is an inverse problem which has the form 

x = M. s + e 

where: x is the data from the telescope i.e. x; is the number of neutrons in the energy bin 
from Ei to E i+ 1; M is the response matrix of the telescope (see Section 4); s is the data as 
it was at the source i.e. what it is that we are actually trying to find; e is the noise on the 
data. It is the noise on the data which makes inverse problems difficult to solve. If there was 
no noise and the problem was just of the form ® 

x = M.s 

then we could simply invert the matrix and the solution would be 

s = M~ 1 .x 
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However, such straightforward inversion amplifies any noise on the data and gives totally 
misleading results. Hence, we must use another method known as statistical regularisation. 
This involves compromising between an exact fit to the data, which is what straightforward 
inversion gives, and a uniform, featureless spectrum. 

The particular statistical regularisation technique used in this project was Maximum En- 
tropyi In order to perform this smoothing technique on the data, the Glasgow University 
Inversion Problem Subroutines (GUIPS) were used. These Fortran subroutines give various 
choices of ‘smoothness’ measure but we concentrated on the Maximum Entropy method. T his 
gave us four possible options, either the Global Maximum Entropy smoothness constraint or 
the Local Maximum Entropy smoothness constraint could be used in conjunction with either 
a chisquared or Bayesian smoothing parameter. The Global Maximum Entropy smoothness 
constraint is of the form 

$(/) = £ f* ~ m i- fi l°g (fi/mi) 

♦ ^ 1=1 

where m* is a ‘prior estimate’ of / which the function is smoothed towards. The Local 
Maximum Entropy constraint is of the form 


*(/) = £(/* - /«+0(log(/i) - log(/ i+ i)) 


i=l 


This modification removes the need for the prior m. To then solve the problem it is necessary 
to find the f which minimises 

X 2 + A$(/) 


where 


* 2 = £ 


*=1 


i 2 


ft - £ Hijfi 
j = 1 


and X is the smoothing parameter. The choice of a Bayesian smoothing parameter puts 
X 2 — X <f> = mcr 2 where m is the number of data points and a 1 is the noise variance on the data. 
Choosing the chisquared smoothing parameter puts x 2 = nw 2 . For this project the GUIPS 
routines were treated as ‘black boxes’ and the actual details of the numerical strategy used to 
optimise x 2 + were not tampered with. 


4 The Response Matrix 

The response matrix Mij gives the probability of a neutron in energy bin j being incorrectly 
counted in energy bin i. The matrix that was used in this project was generated by Monte 
Carlo simulations at the University of California and is shown below. 
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The matrix should be read as follows: 139 neutrons which had a total energy of 13 MeV were 
counted in the first energy bin; 5 were counted in the second bin; 1 was counted in the third 
bin etc. The energy bins which were used are shown in the following table. 
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Bin Energy 

Bin Width 

13 

8.5 - 17.5 

22 

17.5 - 26.5 

32 

26.5 - 37.5 

50 

37.5 - 62.5 

77 

62.5 - 91.5 

100 

91.5 - 108.5 


Approximately 180,000 neutrons at each of the six energies, 13, 22, 32, 50, 77, 100 MeV, 
were ‘injected’ into the detector. Only those events which gave a signal consistent with single 
elastic scattering were retained and the results were binned according to their energies. It was 
assumed that all neutrons arrived at an angle of 15° to the detector axis. This is consistent 
with the large solar flares that occurred in July 1991. 4 The matrix is not a unit diagonal matrix 
because, particularly at low energies, there is a non- negligible chance of a neutron scattering 
m6re than once in the forward detector. This means that the actual energy deposited in the 
forward detector will be underestimated, although the event itself will still resemble an elastic 
scatter event. There is also a possibility that an inelastic scatter event could masquerade as 
an elastic event, giving a spurious result. Also only a small fraction of the incident neutrons 
scatter into the range of solid angles that includes the rearward detector. The response matrix 
is inherently noisy because, as can be seen from the matrix itself, there are only small numbers 
of neutrons in the bins. The matrix was used to generate the fake data required and then the 
data was convolved back through the matrix during the statistical regularisation process in 
order to reconstruct a spectrum so that it could be compared to the original fake spectrum. 


5 Details of the Project 

5.1 Simulation of Data 

The code which this project used had four basic steps. It read in the response matrix and 
other details from a file, generated a fake spectrum, performed statistical regularisation on 
the fake data and then outputted the reconstructed spectrum. Originally, the code generated 
a perfectly flat spectrum using the formula detailed below, 

K v /massy'T--~2e7mass } 

We / 

where K is a constant, mass is the mass of the neutron in MeV and e is the energy of the bin in 
MeV i.e. 13, 22, 32, 50, 77 or 100. However, since this does not correspond to what is actually 
recorded, the first task in this project was to change the code so that the generated spectrum 
better resembled the actual data received from COMPTEL. This meant incorporating the 
decay time of the flare into the calculations for the spectrum so that the spectrum had the 
form of an exponential decay. This is important because the decay times of the flares are of 
the same order of magnitude as the e-folding time of the neutrons, so some of them will have 
decayed before they reach the telescope i.e. observations over a finite time interval monitor 
an energy-dependent interval of the actual flare. The revised formula is 

I tn u\ ^ (Ky/massJ! “ 2e/mass\ ] 
y = exp < ((1 ft) - (1/r)) I — 1 \ 
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where K, mass, e are as above, t is the decay time of the flare and r is the e-folding time of 
the neutron. The flare decay time used during this project was 1000s and the e-folding time 
of the neutron was 930s. The spectrum which was used is shown in Figure 2. This spectrum 
does actually resemble data which is available from the large flares which occurred in July 
1991. 5 

The fake data was also perfectly ‘quiet’ but unfortunately attempts to add substantial 


Figure 2: The Fake Spectrum 



noise directly to the data were unsuccessful. We used a poisson noise routine in C 6 to generate 
random noise which was added to our data points. Because we were adding poisson noise to 
data which was of the order of tens of millions the noise, being the square root of the data, was 
negligible. Reducing the total number of neutrons to increase the noise on the data did not 
help because the Fortran subroutines seemed unable to cope with a total number of neutrons 
below approximately one hundred thousand neutrons, at which stage poisson noise is still 
negligible. Various attempts were made to reach a compromise between non-negligible noise 
and a suitable total number of neutrons or to resolve the difficulties in the Fortran routines 
but this proved to be beyond the scope of this project. 

Instead of adding noise to the data directly, we used the generated data normalised to 
high count rates, to concentrate attention on the uncertainty arising from the poisson noise 
on the response matrix, implicit in the Monte Carlo method. To do this we generated many 
alternative poisson realisations using the response matrix using a shell script. For each of 
these we first convolved the fake data through the response matrix and then reconstructed 
the solution using each of the four possible regularisation methods. However, an unexpected 
difficulty arose: certain realisations of the response matrix resulted in unstable behaviour in 
the GUIPS routine, in which counts in some energy bins were effectively set to zero. See Figure 
3. Investigation of this effect proved to be beyond the scope of this project. We proceeded 
to put the necessary procedures in place by artificially suppressing the response matrix noise 
level by multiplying the numbers ih the noisy response matrix by a factor of 100 and then 
dividing the resulting fake data by the same factor to compensate. This effectively reduced 
the noise level on the matrix by a factor of 10. See Figure 4. A complete investigation of the 
problem was deferred for future work. 
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Figure 3: Using a ‘Noisy’ Response Matrix 



Figure 4: Using a ‘Less Noisy’ Response Matrix 



5.2 Confidence Intervals for the Various Methods 

So that an average spectrum could be plotted, more code was written in C to calculate the 
weighted mean of the total number of neutrons for each energy bin, as well as the length of 
error bar required. The error bars on Figure 5 show the region of 90% probability, i.e. 90 out 
of the 100 realisations gave a value for the total number of neutrons which was between these 
two points for that particular energy bin. Figure 5 must obviously be treated with care due 
to the small counts generated by the Fortran routines in some bins for certain realisations. 


Figure 5: The Mean Spectrum using a ‘Noisy’ Response Matrix 
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5.3 Choice of Regularisation Strategy 

In order to choose which of the methods of statistical regularisation was best, the shell script 
was used to create 100 realisations for each of the four possible choices. We used the noisiest 
response matrix for this so that we could see which of the methods were most sensitive to 
this. Figures 6, 7, 8 and 9 show the mean spectra obtained from these realisations, with 
the error bars showing the region of 90% probability. The original data is also plotted here 
(denoted by a +). From these graphs it can be seen that, while still returning spurious 

Figure 6: Local Maximum Entropy with Chisquared Smoothing Parameter 



Figure 7: Local Maximum Entropy with Bayesian Smoothing Parameter 
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Figure 8: Global Maximum Entropy with Chisquared Smoothing Parameter 
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results in some energy bins for some realisations, the best choice of regularisation method 
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Figure 9: Global Maximum Entropy with Bayesian Smoothing Parameter 



is the Local Maximum Entropy Constraint with a Bayesian choice of smoothing parameter. 
The Global Maximum Entropy smoothing constraint is more sensitive to possible features 
at the edges of the spectra than the Local Maximum Entropy smoothing constraint, which 
means for this particular case that it is more to susceptible to noise in these areas. Choosing 
the chisquared method tends to produce a larger than optimum smoothing parameter, which 
results in the data being over-smoothed. The Bayesian choice selects a smoothing parameter 
which is slightly smaller than the optimum value but in this case this had no adverse effects 
on the final result. 

6 Conclusions 

The main conclusion of this project was establishing the procedures, as detailed above, for 
generating true neutron energy distributions incident on COMPTEL. In particular, we have 
demonstrated that the Local Maximum Entropy ‘smoothness’ constraint, in conjunction with 
a Bayesian choice of smoothing parameter yields the best data fit, and is least sensitive to 
noise on the response matrix. Figure 10 shows the mean spectrum obtained using the Local 
Maximum Entropy constraint with a Bayesian choice of smoothing parameter and the response 
matrix with reduced noise. 

Figure 10: The Mean Spectrum Using the ‘Less Noisy’ Response Matrix 
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7 Further Work 


A further investigation into ways of adding noise to the data itself, while still operating within 
the constraints of the Fortran subroutines, should produce some helpful results on exactly 
how accurate a solution for the inverse problem we can obtain. This would possibly be better 
effected by altering the constraints on the Fortran routines so that they could cope with 
reduced numbers of neutrons, to allow the poisson noise to reach reasonable levels. It will also 
be necessary to investigate why the subroutines give strange results for the first two energy 
bins. Another interesting experiment would be to run the reconstructions with real data from 
COMPTEL to see what sort of features the actual spectrum has. It is important to remember 
that it is not only neutrons which are emitted during flares and it would also be useful to 
collate information gathered on 7 -rays and X-rays with the information available on neutrons 
to' form a more complete picture of what happens in the sun during a flare. However, a lot 
of work still has to be done before an absolutely complete spectrum from a solar flare can be 
obtained. 
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